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Tt ; ABSTRACT 
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. We employ an approximate treatment of dissipative hydrodynamics in three dimensions to study 

the coalescence of binary neutron stars driven by the emission of gravitational waves. The stars are 
^jQ. modeled as compressible ellipsoids obeying a polytropic equation of state; all internal fluid velocities 

^ ' are assumed to be linear functions of the coordinates. The hydrodynamic equations then reduce to a 

set of coupled ordinary differential equations for the evolution of the principal axes of the ellipsoids, the 
internal velocity parameters and the binary orbital parameters. Gravitational radiation reaction and 
viscous dissipation are both incorporated. We set up exact initial binary equilibrium configurations 
and follow the transition from the quasi-static, secular decay of the orbit at large separation to the 
rapid dynamical evolution of the configurations just prior to contact. A hydrodynamical instability 
J> ' resulting from tidal interactions significantly accelerates the coalescence at small separation, leading 

'nI" , to appreciable radial infall velocity and tidal lag angles near contact. This behavior is reflected in the 

' gravitational waveforms and may be observable by gravitational wave detectors under construction. 
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1. INTRODUCTION 



Coalescing neutron star binaries are the primary targets for the detection of gravitational waves 
by the planed LIGO/VIRGO laser-interferometer system (Abramovici et al. 1992). The event rate 
of binary coalescence has been estimated to be ~ 10^ yr~^ Gpc~^ (Narayan, Piran & Shemi 1991, 
Phiniicy 1991). Extracting gravity wave signals from noise requires accurate theoretical waveforms 
in the frequency range 10 — 1000 Hz, corresponding to the last few minutes of the binaries' life and 
orbital separations less than about 700 km (Cutler et al. 1993). 

To leading order, the binary inspiral and the resulting waveform are described by Newtonian 
dynamics of two point masses, together with the lowest-order dissipative effect corresponding to the 
emission of gravitational radiation via the quadrupole formula. One important correction is the effect 
of post-Newtonian terms, which can produce large cumulative orbital phase error (e.g., Cutler et 
al. 1993), accelerate binary coalescence at small separation (Lincohn & Will 1990; Kidder, Will & 
Wiseman 1992, 1993), and cause precession of the orbital plane (Apostolatos et al. 1994). The other 
corrections come from hydrodynamical effects due to the finite size of neutron stars. The analysis of 
Bildsten & Cutler (1992) showed that the binary cannot be synchronized by viscous torque during 
the orbital decay (see also Kochanek 1992). Thus Newtonian spin-orbit coupling has negligible effect 
on the orbital phase at large separation, unless the neutron stars have intrinsic spins near the break- 
up limit (see also Lai, Rasio & Shapiro 1994a, hereafter LRS3). Another aspect of hydrodynamical 
binary interactions concerns resonant tidal excitations of g-modes in the stars, which occur at large 
orbital separation (orbital frequency <, 100 Hz). It was shown (Reisenegger & Goldreich 1994; Lai 
1994), however, that the "resonant tides" and their effects on the orbital decay rate are rather small. 
Therefore, it has now become clear that hydrodynamical interactions are important only during the 
final stage of neutron star binary coalescence, when the orbital separation is within a few stellar 
radii. The importance of tidal effects havs been demonstrated in our previous studies (LRS3; Lai, 
Rasio & Shapiro 1993b, hereafter LRS2), where it was shown that close binary systems containing 
sufficiently incompressible fluid, such as binary neutron stars, are dynamically unstable as a result of 
strong Newtonian tidal interactions. The basic consequence of the instability is the acceleration of 
the coalescence, leading to a significant radial infall velocity at contact, comparable to the free-fall 
velocity. 

At large orbital separation, before the stability limit is reached, the neutron star binary evolves 
quasi-statically along an equilibrium sequence. This phase of the evolution has been studied in detail 
in LRS3, where we have considered the effects of varying neutron star masses, radii, spins and the 
stiffness of the equation of state. We have also shown that stable mass transfer is nearly impossible. 
At small orbital separation, the quasi-equilibrium treatment of LRS3 can provide qualitative features 
of the binary evolution. However, to obtain quantitatively accurate results when the orbital evolution 
takes place on a timescale comparable to the internal hydrodynamic timescale, especially when a 
dynamical instability is encountered, we need to use fully dynamical equations to describe the system. 
A new approximate energy variational formalism based on compressible ellipsoidal figures has been 
developed recently (Lai, Rasio & Shapiro 1994c, hereafter LRS5) to handle binary hydrodynamics. 
The essence of our treatment is to replace the infinite number of degrees of freedom in a fluid binary 
by a small number of variables specifying the essential geometric and kinematic properties of the 
system, such as ellipsoid axes, angular velocity and vorticity. The hydrodynamics is then described 
approximately by a set of ordinary differential equations (DDEs) for the time evolution of these 
variables, in place of the usual hydrodynamical PDEs. In this paper we apply this dynamical theory 
to neutron star binary evolution just prior to the final merging. 

Our dynamical ellipsoid binary model represents the compressible generalization of the Riemann- 
Lebovitz equations for an isolated, incompressible ellipsoid (Lebovitz 1966; see also Chandrasekhar 
1969, hereafter Ch69). Our approximation scheme provides an exact description of Newtonian binary 
stars when (a) the fluid is incompressible (polytropic index n = 0), and (b) the tidal effects beyond 
quadrupole order can be ignored. Hence for neutron stars governed by a moderately stiff equation of 
state (n ^ 0.5) and separated by a few stellar radii, our model is an excellent approximation to the 
true solution. Although our dynamical model is essentially equivalent to the affine model developed 
by Carter & Luminet (1983, 1986) in the context of tidal interactions with a massive black hole, our 
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fornmlation of the' problem is quite difFcrcnt and makes explicit use of global quantities such as the 
total angular momentum and fluid circulation, which are conserved in the absence of dissipation. As 
a result, appreciable simplification in the description of many dynamical processes can be achieved. 
Moreover, we incorporate both viscosity and gravitational radiation reaction as possible dissipation 
mechanisms. 

Complete understanding the final coalescence and merger of the neutron stars requires full 3D 
hydrodynamical simulations. So far, all simulations start from a binary configuration near contact, 
by which point hydrodynamical effects are already important (Oohara & Nakamura 1990; Nakamura 
& Oohara 1991; Shibata, Nakamura & Oohara 1992, 1993; Rasio & Shapiro 1992, 1994; Davies et 
al. 1994). Modeling the distant, pre-contact phase of binary coalescence by 3D numerical codes would 
require prohibitively large amounts of computer resources. Our formulation of binary dynamics in 
terms of ODEs allows us to follow the evolution of the system over a large number of orbits without 
having to worry about excessive computational time or about the growth of numerical errors. Thus 
our dynamical binary model is particularly useful in studying the pre-contact transition phase in which 
the binary evolves from the quasi-static secular regime to the fully dynamical regime. In addition, 
our model can provide reliable initial data for 3D simulations of the binary merger and can serve as 
a check of more sophisticated numerical routines. 

The main purpose of this paper is to present the complete dynamical equations for two stars 
modeled as ellipsoids in binary orbit (§2), including dissipative forces of viscosity and gravitational 
radiation reaction (§3). The equations are used to demonstrate the dynamical instability in the binary 
(§4) and to study a few selected scenarios for neutron star binary coalescence driven by gravitational 
radiation (§5). No attempt is made at a complete survey of parameter space. Instead, we show how 
straightforward it is to vary the individual masses, equation of state and spins of the interacting stars 
so that our method can be easily employed by future investigators to examine other cases of interest. 

We adopt geometrized units and set G = c = 1 throughout the paper. 

2. DYNAMICAL EQUATIONS OF DARWIN-RIEMANN BINARIES 

The dynamical equations for compressible Roche- Riemann binaries, consisting of a finite-size star 
and a point mass, have been derived in LRS5. The readers are referred to that paper for the basic 
definitions of variables, notation and further details. The generalization to compressible Darwin- 
Riemann binaries, consisting of two finite-size stars, is straightforward and is sketched below. 

Consider a binary system containing two stars of mass M and M', each obeying a polytropic 
equation of state. Throughout this paper unprimed quantities refer to the star of mass M and primed 
quantities refer to the star of mass M'. The density and pressure are related by the relation 

P = iiTp^+i/", P' = K'p'^+^Z''' . (2.1) 

Note that for given n and n', the values of K and K' are determined from the equilibrium radii Ro 
and i?Q of the nonrotating stars with the same masses. Thus any realistic equation of state can be 
approximately fit by a polytropic law (cf. LRS3, §4.1). 

In the Darwin-Riemann binary model, both stars have the structure of compressible Riemann-S 
ellipsoids (Ch69, Lai, Rasio & Shapiro 1993a, hereafter LRSl). A general Riemann-S eUipsoid is 
characterized by the angular velocity CI = Cle^ of the ellipsoidal figure (the pattern speed) about a 
principal axis and the internal motion of the fluid with imiform vorticity C = C^z along the same axis 
(in the frame corotating with the figure). We assume that the surfaces of constant density inside the 
star form self-similar ellipsoids. The number of degrees of freedom for star M is then reduced from 
infinity to five: three principal axes ai, a2, as, and two angles (p, ijj, defined such that d(f)/dt = f2 and 
dip/dt = A, where A = —aia2(^/{a\ + a|). Similarly, for star M' , we have a'^, a'2, 03, (j)' and ^' . The 
spins are assumed to be aligned parallel to the orbital angular momentum. The orbit is specified by 
the orbital separation r and the true anomaly 6 (see Fig. 1). Thus, a total of 12 variables completely 
specify the dynamics of the binary system. 

The Lagrangian of the system can be written as the sum of the intrinsic stellar components Lg 
and L'g, and the orbital component Lorh according to 

L = Lg + L'g+ Lorb- (2.2) 
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To derive the stcillar term Lg, let ei, ei and 63 be the basis unit veetors along the instantaneous 
direction of the principal axes of the ellipsoid M, with 63 perpendicular to the orbital plane (the 
'^body frame"). In the inertial frame, the velocity field u of the fluid inside M relative to the center of 
mass can be written as 



u 



ai 



a2 



— A - n X2ei + A + n xie2 



a2 



Oi 



ai 02 , as 

— xiei H X2e2 H 3:363 

ai 02 03 



The kinetic energy of M relative to its center of mass is then given by 

12 1 
T = -/(A^ + n^) - -KnMaia2Ml + — K„M(d^ +al + al), 



(2.3) 



(2.4) 



where I = KnM{al + ai)/^ ^^'^ k„ < 1 is a constant of order unity which depends only on n (see 
Table I in LRSl). 

The internal energy of M is given by 



U = jn^dm = kiKpl/"M. 



(2.5) 



where ki is another constant depending only on n and pc oc 1/(010203) is the central density. The 
self-gravitational potential energy is given by 



W 



3 M-" I 
"5-n R 2i?2' 



with I = Aial + A2al + A^al, 



(2.6) 



where R = (aia2a3)^/^ is the mean radius of the ellipsoid, and dimensionless index symbols Ai are 
defined as in Ch69 (§17). Therefore we have 



Ls=T -U -W. 

Similar expressions for M' can also be derived. The orbital component is clearly 

Lorb = ^fir^ + ^fir^0^ - Wi, 



(2.7) 



(2.E 



where ji = AI M' /{M + M') is the reduced mass. The gravitational interaction energy W, between M 
and M' is given to quadrupole order by 



Wi = - - ^[/ii(3cos2 a - 1) + /22(3sin2 - 1) - 733] 

- ^[^n(3cos' a' - 1) + /^2(3sin2 ^' _ 1) _ j^^j^ 



(2.9) 



where a = 9 — (j), a' = 6 — (f)' , and = K„MafSij /5, and similarly for /'^ . 

Given the Lagrangian, the dynamical equations can then be obtained from the Euler-Lagrange 
equations 

d dL dL 

Mir = 7^' 2.10 

dt dqi aqi 



where {qi} = {oi, (p, ip, a'^, (j)', ip' , r, 9}. For qi = (p, we have 

dJs 3M' 



dt 2r3 

where Jg is the "spin" angular momentum of M 



sin2a(/ii - 122) = AT, 



dL 2 ,^ 

= ^ = ~ -K;„Mai02A, 
oil 



(2.11) 



(2.12) 



and M is the tidal torque on the star. Similarly, for Qi = (j)' , we have 

dJ' 3M 



dt 



2r3 



sin2a'(/(i-/^2)=Ar'. 



For qi = 6, equation (2.13) gives 



dJ, 



orb 



dt 



where Jorh = ^r^O is the orbital angular momentum. Thus the total angular momentum 

J — J s J g ~\~ Jorb^ 

is conserved. Finally, for qi = ip and qi = tp', we obtain 



dC 
'dt 

where C is the fluid circulation in star M 

dL 



= 0, 



d£_ 
'dt 



= 0, 



C = = /A - -KnMaia^Ct, 
oA 5 



(2.13) 



(2.14) 



(2.15) 



(2.16) 



(2.17) 



and similarly for C. Thus the fluid circulations in M and M' are individually conserved. 

The complete dynamical equations can be written in a numerically convenient form as follows: 



til = ai {fP + A^) - 2a2f^A aiAip + 

Qn 



27r 

0,2 = 02(0^ + A^) - 2ainA 02^2/0 + 

Qn 



1 5— (3 cos a - 1), 

TlKn Pc J ai 



27r , _ / 5fci P, 

as = 03^3^+ 

qn \nKn Pc ) as 

-1 r 



M'aa 



0=1—-— 
ai 02 



ft A 



A = / — - — 
ai 02 



02 ai 



n A 



2 — + — di-2 — + — da 



5fci 

niin pc J 0,2 



n A 



M'a2 



(3sin2a- 1), 



3M' 



ai 02 



2 — + — di-2 — + — d2 



ai 02 



O A\ 



ai 02 1 . „ 
1 sm2a 



2r^ \a2 ai 
3M 



■ sin 2a 



0,2 o\ J 



r = rO 



M + M' 3k„ M + M' 
r2 10 ^ 



[ai(3 cos^ a - 1) + 05(3 sin^ a - 1) - Og] 



[af (3 cos^ a' - 1) + 4^(3 sin^ a' - 1) - a'i] 



3k'„ M + M' 

To ^ 

2re 3k„ M + M' 2 2x • o M + M' ,2 ,2^ . ^ / 
-3 {a\ - 4) sm 2a - ^ — (af - a'i) sm 2a', 
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(2.18) 

(2.19) 
(2.20) 

(2.21) 
(2.22) 

(2.23) 
(2.24) 



where 9 = florb^ 4> — and g„ = k„(1 — n/S). The equations for dj, Cl' and A' can be similarly written 
down by switching unprimed quantities with primed quantities in equations (2.18)-(2.24). Also, the 
pressure term {bk\Pc) / {nKnpc) can be conveniently expressed in terms of Ro, M and other dynamical 
variables as 



5fci Pc 



M {Ro 



3/ra 



(n ^ 0), 



ni^n pc QnRo \ R 

(see LRS5, where the equivalent expressions for the limiting case of n = can also be found). 



(2.25) 



3. GRAVITATIONAL RADIATION REACTION AND 
VISCOUS DISSIPATION 

Dissipation modifies Euler-Lagrange equations according to 

d_dL^ _dL_ ^ 
dt dqi dqi ^' ' 



(3.1) 



where Tq^ is the "generahzed" force associated with variable gj, and it is defined via the non- 
conservative dissipation rate W (Rayleigh's dissipation function, cf. Goldstein 1980) by W = J^qSi- 

3.1 Gravitational Radiation Reaction 

The main driving force for neutron star binary coalescence is gravitational radiation reaction. 
Consider the orbital coordinates {e^}, centered at the center-of-mass (CM) of the system with ej 
along the line joining MM' , 83 pcrpcndcular to the orbital plane, and 82 perpendicular to ej and 
eg (see Fig. 1). This coordinate basis rotates with angular velocity ^orb with respect to an inertial 
coordinate system. In the weak-field, slow motion regime of general relativity, the gravitational 
radiation emission induces a back-reaction potential ^ react, which can be written as (Misner, Thorne 
& Wheeler 1973): 

<^react = li^^iX-j. (3.2) 



where li? is the fifth derivative of the reduced quadrupolc moment tensor of the system projected 
onto the orbital frame. The contribution from M to W is given by 

Wm = - / V • V^reactdm 

Jm 



:K„Mj ^Ijj' cos^ a + 1^1' sin^ a — sin 2a^ aidi 

, /t(5) . 2 , t(5) 2 , t(5) ■ o \ - I • 

-|- ( ijj sm a + cos a + sm2a 1 0202 + igg 0303 



(3.3) 



cos 2a + (1^°/ - 1^^^)^ sin 2a ) (af - a^)n 



r(5) t(5) 



/ 



+ 0r'^ ft h 

^ ^12 ' cm^'-orb 



where v = u -|- Uorb is the fluid velocity in M, rem is the distance from the CM of the system to the 
CM of M. A similar expression can be written down for M'. Therefore 

W =Wm + Wm' 



5 V 5 



^1 cos^ a + 1^1^ sin^ a — sin 2a^ a\a\ 

((5) 2 (5) 2 (5) \ (5) 

Ijj sin a + 1^2 cos a + sin 2a 1 020(2 + 133 03013 

+ (i-ll cos 2a + (Ipj^ - 4f )i sin2a^ {a\ - ai)0 

(4f cos2 a' + 4f sin2 a' - sin 2a') a'^ai 



I { In'M' 



5 V5 " 



(3.4) 



+ (4t^ sin^ a' + 4f cos^ a' + if^^ sin 2a') a'^a'^ + 4g ^303 



ig cos 2a' + - l|0^ sin2a' ) {a'^ - a'i)^' 



-ln(lf^rr + lf^r^norb). 



The dissipativc forces due to gravitational radiation are then given by Tq^ = dW/dcji. Clearly, 
~ ~ 0, thus gravitational radiation reaction conserves C and C. With the inclusion of 

gravitational radiation reaction forces, the dynamical equations (2.18)-(2.24) for binaries are modified 

according to: 



ai 



a2 



03 



n = 



{•••} 
{•••} 
{•••} 

02 _ 

a-L 



Ij cos^ a 



- ii^l^ sin^ a 



f sin 2a 



1(5) • 2 , 1(5) 2 , 1(5) . o 
Ijj sm a + Igg cos a + Ijg sin 2a 



ai, 



a2, 



^1(5) 
5^53 «3, 



Oi 
02 



{•••} + i ( lis cos 2a + 9 - l^J ) sin 2a 



(5) T(5)^ 



Oi _^ 02 

a2 ai 



02 
ai 



A= 



{•••} 



Oi 

02 

2t(5) 



{•••} + - 



ll|^cos2a+^(I,, 



(5) 



■ I^l^ ) sin 2a 



(3.5) 
(3.6) 
(3.7) 

(3.8) 
(3.9) 

(3.10) 
(3.11) 



where {• • •} denote the nondissipative terms that already exist in equations (2.18) (2.24) (This nota- 
tion will be used throughout the paper). The expressions for aj, 17' and A' are similar. 

When the timescalc for both orbital and internal quantities to change is much longer than the 
rotation period, e.g., \dai/dt\ << \nai\, simple expressions for I^^"* can be derived. To order dai/dt, 
the nontrivial components of i^-) are 



(3.12) 



if^ = - = 16f2^(7ii - /22) sin2Q: + - I'^^) sin2a' 

+ AQ^l^^,[2norb^lrf + 2tli^ir^] + 40f7-^[fi(/ii - 122) + 2f2(/ii - 722)] cos 2a 
+ 40Q'3[17'(/(i - 122) + 20' (/(i - 7^2)] cos 2a', 

Ip-) =lg) = + 1617^(711 - 722) cos 2a + - I'^^) cos 2a' 

- m9!^[n{hi - 722) + 20(7ii - 722)] sin2a 

- 40fi'3[17'(7;i - 7^2) + 2f2'(7;i - 7^2)] sin 2a'. 

These generalize the expressions derived in Appendix A of LRS5 for Roche-Riemann binaries. 

3.2 Viscous Dissipation 

Viscous dissipation forces can be easily incorporated into our dynamical equations. Since the 

motion of the center of mass of the star is not affected by viscous dissipation (which depends only on 
the shear stresses inside the star), the viscous forces for an isolated star as derived in LRS5 (§4.1) can 
be directly applied to binaries. The dissipation rate due to shear viscosity is W = Wm + Wm' , with 



4 

Wm = - ^vM 



+ 



02 
02 



+ 



aia2 



(3.13) 



where v is the mass-averaged shear kinematic viscosity v = j vdm/M. The expression for Wm' can 
be written down similarly. Since W is quadratic in 5,, the dissipative forces are given by = 



(l/2)dW / dqi. Clearly, in the presence of viscosity, the fluid circulations C and C arc not conserved. 
However, viscous forces do not affect angular momentum, i.e., !Fe = = J-"^' = 0. The dynamical 
equations (2.18)-(2.24) are modified according to 



, - 10 _ / 2di d2 

oi ={•••} — - — 

3k„ V ai a2 



a2 = {---} 
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as = 



•••}- 




a2 




ai 


02 ) 


0.2 




ai 


02 ) 



202 
02 

/2a3 

V 03 



Ol 
Oi 
Ol 
Ol 



1 

Ol' 

1 

02' 
1 

03' 



-1 r 



10 c 



2 2 

ofo^ 



A= 



5 of — o: 
■} + —i>^ ■■ 

Kn ai02 



(3.14) 
(3.15) 
(3.16) 
(3.17) 
(3.18) 



and the expressions for a^, O' and A' are similar, while the expressions for f and 6 are unchanged. 

4. DYNAMICAL INSTABILITY IN DARWIN-RIEMANN BINARIES 

The dynamical instability in a binary system results from tidal interactions between the two stars. 
The height h of the tidal bulge raised on star M by its companion M' is of order h ~ R{M' /M){R/r)^ . 
This tidal deformation makes the gravitational interaction between M and M' more attractive, and 
the tidal potential energy is 

M'Q M''^R^ 

- (4.1) 



Wtide ~ 



(see eq. [2.9]), where Q ^ KnMRh is the quadrupole moment of M. Thus at sufficiently small binary 
separation, assuming k„ is not too small (i.e., the star is not too compressible), the binary interaction 
potential energy ~ {—MM'/r + Wude) becomes steeper than the point-mass contribution —MM'/r. 
This is the cause of the dynamical instability in the binary, as is common to all interaction potentials 
that are sufficiently steeper than 1/r (cf. Goldstein 1980, §3.6). 

To determine the dynamical stability limit, one only needs to construct a sequence of equilibrium 
models with constant C and C . The onset of dynamical instability corresponds to the turning point 
in the energy curve along the sequence. This procedure is described and applied in LRSl and LRS4, 
where numerical values of the stability limits for various binary models have been tabulated. 

Using our dynamical equations with no radiation redaction or viscosity, we can now study how 
the instability develops in time. In Figure 2 we show an example of the time evolution of an unstable 
system with n = n' = 0.5, K = K', p ^ M/M' = 1/2 and A = A' = (corotation). The dynamical 
stability limit is at f = r/{ai + a'l) = 1.174 (cf. LRS4, Table 2). The dynamical equations are 
integrated numerically using a standard fifth-order Runge-Kutta scheme with adaptive stepsize (Press 
et al 1992). At t = 0, an unstable equilibrium solution is constructed for f = r/(oi +a'i) = 1.17. This 
equilibrium solution is then perturbed by setting f — — 10^^(M/i?o)^^^- For comparison, the results 
of an integration for a stable binary with f = 1.18, and with the same applied perturbation are also 
shown. We see clearly that as the dynamical instability develops, oi increases while r decreases, and 
this is accompanied by the significant development of tidal lag angle in the two stars (a, a' > 0) and 
de-synchronization (A, A > 0). Of course, the precise evolution of an unstable binary depends on 
how the initial configuration is perturbed. 

The development of tidal lag in the absence of fluid dissipation can be qualitatively understood 
as follows (cf. Lai 1994). For star M in the binary system, the tidal potential cx 1/r^ due to the 
companion M' acts like an external perturbing force, with a driving frequency ~ Ail = VLorb — 
where ^orb and are the orbital and spin angular frequencies. The star has an intrinsic dynamical 
frequency of order ~ (M/i?^)^/^. Schematically, the equation governing the tidal distortion ^ can 
be written as 

■■ " ^ -"'^"^ (4.2) 
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If r and Af2 were constant in time, then the stationary tide would be; in phase with the driving force. 
However, when r and/or florb change(s) during binary evolution, we have (assuming lOo >> AJ7) 

- ladyn) 



where the lag angle is given by 



An , , 

Oidyn ^ (4.4) 



and td ~ \r/f \ ~ \norbl^orb\ is the orbital decay timescale. Thus a dynamical tidal lag arises naturally 
even without fluid dissipation, and it is due to the finite time necessary for the star to adjust its 
structure to the rapidly changing tidal potential. Since the orbit decays rapidly as a result of dynamical 
instability, the binary de-synchronizes, and thus adyn becomes large at small r. 



5. NEUTRON STAR BINARY COALESCENCE 

The main parameters that enter the evolution equations of §2 and 3 are the mass ratio p = M/M' , 
and the ratios Rg/M and R'^/M' . The ratio Rq/M is determined from the nuclear equation of state 
(EOS). For the canonical neutron star mass M = IAMq, all EOS's tabulated in Arnett & Bowers 
(1977) give Ro/M in the range of 4-8. Recent microscopic nuclear calculations indicate that the 
neutron star radius Rq is ~ 10 km, almost independent of the mass for M in the range of O.SM© to 
1.5M0 (Wiringa, Fiks & Fabrocini 1988). Thus we will use Rq/M ~ 5 as a representative value and 
assume Ro = R'o- 

A polytrope is only an approximate parametrization for a real EOS. We can determine the 
approximate polytropic index n that mimics the structure of a real neutron star by using the tabulated 
values of the moment of inertia of neutron stars (see LRS3, §4.1). Typically we find n ~ 0.5 for 
M ~ 1.4Mq for the EOS of Wiringa, Fiks & Fabrocini (1988). 

Other parameters needed for the calculations are the spins of the two neutron stars fig and 
when the binary separation is large. For a uniformly rotating neutron star, the spin rate is limited 
by mass-shedding, = 0^/(M/i?3)i/2 ^ o.6 (Friedman, Ipser & Parker 1986, Cook, Shapiro & 
Teukolsky 1992). At large orbital separation, we have ai — > 02 and O — > Q.orb (M-|-M')^/^/r^/^ << 
Og, thus Js —IK and C ^ IK (see eqs. [2.12], [2.17]). The fiuid circulation inside the star, which is 
conserved during the binary evolution in the absence of viscosity, is therefore given by C = — 7J7s. Here 
we identify ^Ig = — A(r = 00) as the spin angular velocity at large r (for an axisymmetric spheroid, 
uniform spin and vorticity are indistinguishable; cf. eq. [2.3]). Note that when fig is positive (the spin 
is in the same direction as the orbital angular momentum), C is negative. For simplicity we assume 
fl'g = = C. 

In Figure 3 we show two examples of the pre-contact evolution of binary neutron stars. Both 
calculations start at t = with r = 5Ro, and terminate when the surfaces of the stars contact. The 
coalescence is driven by gravitational radiation reaction, and we ignore viscosity for now. In the first 
example, the stars have equal masses, with Ro/M = R'^/M' = -5, n = n' = 0.5, and both have zero 
spin fig = fig = 0. The initial configuration is obtained by solving the Darwin-Riemann equilibrium 
equations (see LRS4) with r/(ai + a[) = 2.461. In the second example, we have n = n' = 0.5, but 
now M/M' = 1.2 so that Ro/M = 5 and R'o/M' = 6. Also we set Og/(M/i?2)i/2 = o.4 (near the 
maximum possible value) and fi'^ = 0, corresponding to C/{M^RoY/'^ = -0.1204 and C = 0. The 
initial state has r/{a\ -\-a'i) = 2.385 (and r/Ro = 5), and the nondimensional vorticity parameters are 
fa = —{a{ + a\)K/{aia2n,) = 3.818, — —2. Also shown in Figure 3 are the energy curves of the 
constant— C equilibrium sequences. Initially, the binary closely follows the equilibrium sequence. As 
the dynamical instability develops, the radial velocity increases considerably, and thereafter the two 
stars merge hydrodynamically in just a few orbits. The terminal radial velocity at contact typically 
reaches 10% of the free-fall velocity. This qualitative behavior has already been observed in the 
simplified calculations we presented in LRS2 and LRS3, where the stars were constrained to retain 
their equilibrium shapes throughout the evolution. 
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From Fig. 3, wc also sec; that the tidal lag angle a or a' increases with decreasing r, attaining a 
large value of about 10° at binary contact. Let us consider the lag angle a in star M. In the absence 
of fluid viscosity, there are two contributions to this lag angle, a = adyn + ctgr: one is the dynamical 
lag ctdyn discussed in §4, which arises from the rapid orbital decay (especially when the dynamical 
instability develops) ; the other is the gravitational radiation dissipation lag agr analogous to the usual 
viscous lag (see later). During the secular evolution phase (before the dynamical instability develops), 
the orbital decay timescale is 

where Mt = M + M'. Prom equation (4.4), the dynamical lag is given by 

RlAn / 64 MM' MA MRIM'mI'"^ , . 

where in the second equality we have specialized in the case when fig ~ and Af2 ~ Qorfc — 
{Mt/r^y^"^. The gravitational radiation dissipation lag arises because gravitational radiation reac- 
tion directly exerts a torque on the star. Prom equation (3.4), this torque Afgr = = dW/dcf) is 
given by 

0/1 \ r 1 1 



2 fl 



ig^ cos 2a + (4i^ - 4f ) \ sin 2a 



^ (^Kr,M ] (l6Q^^bM?-^ cos 2a) {a\ - al), 



(5.3) 



where in the second equality we have used equation (3.12). Since the fluid circulation C in the star is 
conserved in the absence of viscosity, and since \ Js\ — \C\ to the leading order (cf. eqs. [2.12], [2.17]), 
there must be a small misalignment agr of the tidal bulge so that Mgr can be balanced by the tidal 
torque N (cf. eq. [2.11]). Requiring dJg/dt = Af + Mgr = 0, we obtain 

32 MM.^/^ 

which is independent of the radius and spin of the star. This result was derived previously in LRS5 
(§9.2) for Roche-Riemann binaries. Although agr is larger than adyn for large binary separation 
r ^ 2Ro{M' /My^^, the radiation dissipation lag is always small {agr ^ 0.01). At smaller orbital 
separation, the dynamical lag dominates, Moreover, as the dynamical instability develops, the orbital 
decay time becomes comparable to the orbital period, td ~ torb = ^/^orb, and the thus the lag angle 
near binary contact is 

Por M ~ M' and r ~ 2.6i?o (contact), we have a ~ 0.2, in agreement with the numerical results of 
Pig. 3. Equation (5.5) implies that a is smaller for the spinning, more massive star, also in agreement 
with Fig. 3. 

The gravitational wave forms can also be calculated using the standard quadrupole formula, 
giving 

h+ = [Atr2^]2^^cos26' + (/ii -/22)f^^cos2(/)+(/(i - cos2(t)'] (1 +00329), 

? (5.6) 
h^=-- [^lr'^Q.l^^^sm2e + {hi - hi)^^ sm2(j) + {I[^ - I'^^)Q.''^ sm2(l)'] cos 6, 

where D is the source distance, and 6 specifles the angle between the direction of wave propagation 
and the 2:-axis. We have neglected higher-order terms which are smaller by a factor of order \r/rQ.orb\- 
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Equation (5.6) generalizes the expressions derived in LRS3 to the case when cf) ^ 6 and (j)' ^ 6 
(cf. Fig. 1). Neglecting the small tidal correction, the wave frequency is / ~ (M + M'Y^'^r~^^'^ /it = 
1123 Mf4i(i?o/5Af)-''/2(j./3ij^)-3/2 (for M = M'). Figure 4 depicts the wave amplitude seen 
by an observer along the rotation axis (O — 0) for one of the cases (O^ = 17'^ = 0) shown in Figure 
3. Comparing to the result for two-point masses, we see that a large phase error (2 — 3 cycles) 
develops during the last ~ 10 wave cycles as the result of the accelerated coalescence induced by the 
dynamical instability. Such signature should be detectable by specially configured ("dual recycled") 
interferometers that operate over adjustable, narrow bands around 1000 Hz (e.g.. Strain & Meers 
1991). 

The dynamical effects of viscous dissipation can also be incorporated in the calculations. Dimen- 
sionally, the maximum possible value of viscosity is Pmax ~ {MRoY^^. However, the viscosity due 
to electron-electron scattering (Flowers & Itoh 1979), which is the dominant source of microscopic 
viscosity since neutrons and protons are likely to be in superfluid states in cold coalescing neutron 
stars, is many orders of magnitude smaller than Pmax- Nevertheless, to illustrate the qualitative 
effects, we consider an extreme example, adopting u = O.bRoiM/RoY^'^. In Figure 5, we compare 
neutron star binary evolution with and without viscosity. We see that viscous dissipation tends to 
synchronize the binary, i.e., to reduce A as compared to the inviscid cases. Such synchronization is 
the result of viscous tidal lag a^g ~ ACl/{ui'^tvis) (compare with eq. [4.4]), where tms ~ Ro/^ is 
the viscous timescale (e.g., Zahn 1977; also see eq. [8.4] in LRS5). Also viscous dissipation tends to 
accelerate the coalescence. This is because orbital angular momentum is transferred to the stellar 
spin via viscous torque. However, as can be seen In Fig. 5(c), even with such an extreme value of 
viscosity, synchronization can never be achieved. In fact the binary always becomes more and more 
asynchronized as r decreases, even if it is corotating at large separation. This conclusion agrees with 
that of Bildsten & Cutler (1992). 

Finally, it should be noted that our treatment so far has ignored post-Newtonian (PN) effects 
other than the lowest-order dissipative effect corresponding to the emission of gravitational radiation. 
However, for the typical value of Ro/M = 5, other PN effects are important and can alter the 
orbits considerably. In particular, even in the case of two point masses, these PN effects can by 
themselves make a circular orbit become unstable when the separation is smaller than some critical 
value ("inner-most stable orbit") rcR- Kidder, Will & Wiseman (1992) have recently obtained rcR — 
6(M + M') -I- 4/i. The PN effects lead to a plunge orbit for r < ran., with significant infall radial 
velocity. For two point masses with M = M', ran — 14M. Since the Newtonian hydrodynamical 
stability limit is at ~ iRo typically, tgr and are comparable for Ro/M = 5. Clearly, without 
including PN terms, our final coalescence trajectory can only be approximate. However, it is clear 
from our discussion that the Newtonian hydrodynamic effects discussed in this paper are at least as 
important as relativistic corrections to the orbital motion for the final phase of neutron star binary 
coalescence. When the two effects are both properly incorporated, the final coalescence is likely to be 
even faster, and assume a significant "head-on" character. 

6. CONCLUSIONS 

We have presented a simplified hydrodynamical treatment of close binary systems based on com- 
pressible ellipsoids obeying a polytropic equation of state. We employed this treatment to demonstrate 
the development of dynamical instability during the final phase of neutron star binary coalescence 
prior to contact. Such instability is accompanied by a significant radial infall velocity and dynamical 
tidal lag angles of order 10°. We have also shown that the neutron star becomes more asynchronized as 
the binary orbit shrinks (see also Bildsten & Cutler 1992). Full hydrodynamical simulations of binary 
merger based on non-corotating models with large initial plunging velocity components (Shibata el 
al. 1992, 1993) show significantly different results from simulations based on corotating binary models 
without large infall velocities (Nakamura & Oohara 1991, Rasio & Shapiro 1992, 1994). However, 
these simulations should be considered only preliminary since the calculations start near contact, by 
which time dynamical effects are already important and the initial conditions in these simulations are 
very approximate. 

Our dynamical binary models provide a new tool to study the pre-contact evolution of binary 
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neutron stars. The simplicity of replacing the full hydrodynamical equations with ODEs allows us 
to sample a large number of physical parameters (e.g., masses, spins, equations of states, etc.) and 
incorporate dissipative effects (gravitational radiation and viscosity) easily. Our Newtonian dynamical 
binary equations, when coupled with the post-Newtonian equations for point-mass binaries (e.g., 
Kidder et al. 1993), should give a reasonable description of the terminal phase of the binary prior to 
merger. 

The numerical codes implementing the equations presented in this paper can be obtained from 
the authors upon request. 
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Figure Captions 



FIG. 1. — The coordinate system used for Darwin-Riemann binaries. 

FIG. 2. — The development of the dynamical instability in a binary system with p — M/M' — 1/2, 
n = n' = 0.5, K = K' . The initial configurations are in equilibrium, and corotating. The left panels 
show the evolution of a dynamically unstable binary with f = r/{ai + a'l) = 1.17, while the right 
panels show that of a stable binary with f = 1.18. The solid lines correspond to star M, the dashed 
lines correspond to M' . Each line is labeled by the value of Hn Oj. 

FIG. 3. — The evolution of coalescing binaries driven by gravitational radiation. The system energy 

E, tidal angle a, radial velocity Vr = r and time t arc shown as a function of binary separation r. 
Here n = n' = 0.5, Ro/M = 5, Rq = R'^, and the calculations start {t = 0) at r/Ro = 5 and terminate 
when the surfaces of the stars contact. The left panels show the case with M = M' and = = 0; 
the right panels show the case with M/M' = 1.2, ns/iM/RlY^"^ = 0.4 and n'^ = 0. The dotted curve 
in the E{r) diagram is the equilibrium energy of a constant-C binary sequence. The turning point 
marks the onset of dynamical instability. The dashed line in the right panel corresponds to a', while 
the solid line corresponds to a. 

FIG. 4. — The waveform from neutron star binary coalescence shown in Fig. 3 (fig = OJ^ = 0). The 
dark solid line corresponds to zero viscosity, the light solid line assumes P = 0.5(Mi?o)^^^. The dotted 

line is the result for two point masses. 

FIG. 5. — Radial velocity Vr, number of orbital cycles Ngrb and vorticity parameter A (measuring 
the degree of non-synchronization) of a coalescing neutron star binary. Here M = M', Ro/M = 5 
and n = 0.5. The solid lines arc for C = C = (irrotation) with zero viscosity, the dashed lines 
for C = C = at r/Ro = 5 and with v = Q.^{MRofl'^] the long-dashed lines are for A = A' = 
(corotation) at r/Ro = 5 with no viscosity, and the dotted-dashed lines are for i? = 0.5(Miio)^/^. The 
dotted lines in (a)-(b) are the results for two point masses. 
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